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ABSTRACT 

We present the results of five years (2005-2009) of MAGIC observations of the BL Lac object 
PG 1553+113 at very high energies (VHEs, E > 100 GeV). Power law fits of the individual years are 
compatible with a steady mean photon index T — 4.27 ± 0.14. In the last three years of data, the 
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flux level above 150 GeV shows a clear variability (probability of constant flux < 0.001%). The flux 
variations are modest, lying in the range from 4% to 11% of the Crab Nebula flux. Simultaneous 
optical data also show only modest variability that seems to be correlated with VHE gamma ray 
variability. We also performed a temporal analysis of (all available) simultaneous Fermi/LAT data 
of PG 1553+113 above lGeV, which reveals hints of variability in the 2008-2009 sample. Finally, 
we present a combination of the mean spectrum measured at very high energies with archival data 
available for other wavelengths. The mean spectral energy distribution can be modeled with a one- 
zone Synchrotron Self Compton (SSC) model, which gives the main physical parameters governing 
the VHE emission in the blazar jet. 

Subject headings: radiation mechanisms: non-thermal, gamma-rays: observations, BL Lacertae ob- 
jects: individual (PG 1553+113) 



1. INTRODUCTION 

The majority of extragalactic 7 ray sources, both at 
GeV energies and above 100 GeV, are blazars, radio- 
loud active galactic nuclei with a relativistic jet point- 
ing towards the Earth. Their emission is dominated by 
the non-thermal continuum pro duced within the jet an d 
boosted by relativistic effects (|Urrv fc Padovanil fl995 % ). 
The spectral energy distribution (SED) displays two 
broad peaks, widely interpreted as due to synchrotron, 
low frequency peak, and inverse Compton, high fre- 
quency peak, mechanism (although the high energy peak 
could al so be the result of hadronic processes, as pro- 
posed in iMannheiml fl993l ) . Among blazars, BL Lac ob- 
jects are characterized by extremely weak emission lines 
in their optical spectra, which often makes a measure- 
ment of their redshift difficult. The large majority of ex- 
tragalactic sources detected above 100 GeV are BL Lac 
objects, in which the peak of the synchrotron bump is lo- 
cated in the UV- X-ray bands and the high energy peak 
around 100 GeV (these sources are often called high fre- 
quency peaked BL La cs, HBLs). PG 15 53+113 is a 
BL Lac discovered by iGreen et al.l (|1986l ). The large 
X-ray to radio flux ratio makes this source a typical 
HBL. Indeed, its synchrotron peak is located between 
the UV and X-ray bands. Its optical spectrum is fea- 
tureless, preventing the direct determination of the red- 
shift. Indirect methods based on the non detection of 
the characteristic lines and of the hos t galaxy provide 
lower limits, ranging from 0.09 to 0.78 dSbarufatti et al.l 
[200l ISbarufatti. Treves, fc Falom o 2005). The most re- 
cent estimate, bas ed on the Ly alpha for est method, gives 
a z ~ 0.40-0.45 (|Danforth et al]l2mot) . 

PG 1553+133 has be en discovered as a VH E 7 
ray emitt er by H.E.S . S. (lAharonian et al.l l2006a|) and 
MAGIC (jAlbert et al.l |2007a|), with a flux of approxi- 
mately 2% of that of the Crab Nebula above 200 GeV. 
The spectrum appears extremely soft (photon index 
r ~4), as expected by the absorption of VHE photons 
through interaction with the Extragalactic Background 
Light (E BL) if the source is located at rela tively large 
redshift (jStecker. de Jager. fc Salamonl Il992) . The ab- 
sorption process, in fact, is a function of the energy of 
the photon and of the distance it has traveled. Spectra 
with indices T ~ 4 have been observed in blazars located 
at redshifts above 0.2. 

VHE 7 ray observations have been used as alter- 
native method t o cons train the distance of blazars. 
lAharonian et a l. (2006b) proposed a way to set an up- 
per limit on the distance of blazars based on the as- 
sumption that the VHE intrinsic spectrum, obtained by 



correcting the observed spectrum for the extragalactic 
background light absorption, cannot be harder than a 
fixed value given by theory. The technique, applied to 
PG 15 53+113, lead to an up per limit of z < 0.74. Re- 
cently, Prandi ni et al.l (|2010l ) extended this method us- 
ing the spectrum measured at lower energies as limiting 
slope for the original spectrum, obtaining z < 0.66 for 
PG 1553+113, at 2 a level. Other approaches require the 
absence of a pile up at h igh energies. With this method, 
IMazin fc Goebell (|2007D get z < 0.42. Hence, in the 
case of PG 1553+113, the upper limits obtained with 
these methods are in the range of the limits set by op- 
tical measurements. In this work we adopt the redshift 
z = 0.40. Such a large redshift is also supported by the 
absence of significant points at energies above 700 GeV 
in the spectrum of the source. 

At MeV-GeV energies, PG 1553+113 was not detected 
by EGRET, but it is well visible by the Large Area Tele- 
scope (LAT) on-board Fermi, being detected with a sig- 
nificance above 10 a already in the first three months 
of observations (jAbdo et al.ll2009al) . lAbdo et al.l (|2010b[ ) 
show that the Fermi/LAT spectrum is surprisingly con- 
stant both in normalization and slope over ~ 200 days. 
Interestingly, a stability of the spectrum was also sug- 
gested by the H.E.S.S. and MAGIC observations, show- 
ing rather marginal variability during 2005 and 2006 ob- 
servations. The stability of the VHE 7 ray emission is 
in contrast to the behavior commonly observed in other 
TeV emitting BL Lacs, showing rather pronounced vari- 
ations at all timescales. 

After its discovery, PG 1553+113 was regularly ob- 
served by MAGIC. In this paper, we present the analy- 
sis of the new data taken from 2007 to 2009, combined 
with previous observations. The differential and integral 
fluxes are analyzed, in comparison with partially simulta- 
neous measurements at other wavelengths, and the sta- 
bility of the spectrum over this long period is studied. 
Finally, we combine all the data available and model the 
SED with a one-zone Synchrotron Self Compton (SSC) 
model, constraining the main physical parameters that 
govern the VHE emission in the blazar jet. 

2. MAGIC OBSERVATIONS AND DATA ANALYSIS 

Since autumn 2009, MAGIC (jCortina et al.l I2009D is 
a stereo system composed by two Imaging Atmospheric 
Cherenkov Telescopes (IACTs) located on La Palma, 
Canary Islands, Spain (28.75°N, 17.89°W, 2240m asl). 
In this paper, we present only data collected before 
the stereo upgrade, w ith a single telescope, MAGIC I 
(|Baixeras et al.l I2004D . hereafter called MAGIC. The 
parabolic-shaped reflector, with a total mirror area of 
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236 m 2 , allows MAGIC to collect the Cherenkov light 
and focus it onto a multi-pixel camera, composed of 
577 photo- multipliers. MAGIC camera and trigger are 
designed to record data not only during dark nights, 
but also under moderate light conditions (i.e. moder- 
ate moon, twilight). Due to its comparatively low trig- 
ger energy threshold of ~ 50 GeV, MAGIC is well suited 
to perform multiwavelength observations together with 
instruments operating in the GeV range. 



TABLE 1 
PG 1553+113 FINAL DATA SET 



Cycle 


Date 


Elf. Time 


Zd 


Rate 


DC 






[min] 


[° 


] 


[Hz] 


[M] 




23/03/2007 


58 


19 - 


29 


164 


dark night 




19/04/2007 


32 


22 - 


28 


155 


dark night 




20/04/2007 


150 


17 - 


29 


163 


dark night 


III 


21/04/2007 


115 


17 - 


24 


155 


dark night 




22/04/2007 


101 


17 - 


23 


162 


dark night 




23/04/2007 


143 


17 - 


27 


161 


dark night 




24/04/2007 


92 


17 - 


23 


160 


dark night 




17/03/2008 


58 


17 - 


19 


150 


dark night 




18/03/2008 


26 


18 - 


19 


150 


dark night 




01/04/2008 


43 


20 - 


26 


167 


dark night 




05/04/2008 


109 


17 - 


31 


167 


dark night 


IV 


13/04/2008 


97 


17 - 


22 


117 


dark night 




29/04/2008 


44 


27 - 


36 


151 


dark night 




03/05/2008 


24 


26 - 


31 


146 


dark night 




04/05/2008 


40 


28 - 


36 


155 


dark night 




05/05/2008 


38 


26 - 


33 


150 


dark night 




07/05/2008 


40 


28 - 


36 


153 


dark night 




16/04/2009 


93 


17 - 


27 


133 


2.8 - 3.7 




17/04/2009 


103 


17 - 


28 


151 


1.6 - 2.4 


V 


18/04/2009 


126 


17 - 


28 


168 


0.7 - 1.7 




20/04/2009 


73 


19 - 


35 


171 


0.9 - 1.4 




21/04/2009 


57 


23 - 


34 


177 


0.8 - 1.0 




15/06/2009 


57 


24 - 


35 


125 


0.8 - 2.1 



a PG 1553+113 data set from 2007 to 2009 used in this study. 
From left to right: MAGIC Cycle of observation, first column, and 
corresponding dates in dd/mm/yy, second column; effective time 
of observation in minutes and zenith angle range in degrees, third 
and fourth column. In the last two columns, the rate of the events 
after the image cleaning, in Hz, and the mean DC current in the 
camera, in unit of fiA, are shown. The night is considered as dark 
night, if the DC current while observing an extragalactic object is 
less than indicativcly 1.2 ^A. 

The total Field of View (FoV) of the MAGIC camera 
is 3.5°, and the effective collection area is of the order 
of 10 5 m 2 at 200 GeV for a source close to zenith. The 
incident light pulses are converted into analog signals, 
transmitted via optical fibers and digitised by 2 GHz fast 
analog to digital converters (FADCs). 

PG 1553+113 was observed with the MAGIC tele- 
scope for nearly 19 hours in 2005 and 2006 (jAlbert et al.l 
l2007al) ; it was also the subject of a multiwavelength 
campaign carried out in J uly 2006 with opt ical, X-ray 
and TeV 7 ray telescopes ([Albert et al.l 120091) . Here, we 
present the results of follow-up observations, performed 
for 14 hours in March-April 2007, for nearly 26 hours 
in March-May 2008 , some of those sim ultaneously with 
other instruments (jAleksic et al.l l2010f) , and for about 
24 hours in March- July 2009, which were partly taken in 
moderate light conditions (moon light). Unfortunately, 
both 2008 and 2009 observations were severely affected 



TABLE 2 
PG 1553+113 SIGNAL 



Year 


Time 


Opt. PSF 


Energy Th. 


Excesses 


Signif. 




m 


[mm] 


[GeV] 




M 


2007 


11.5 


13 


80 


1400 ± 242 


5.8 


2008 


8.7 


13 


150 


542 ± 69 


8.1 


2009 


8.5 


14.8 


160 


212 ± 52 


4.2 



a PG 1553+113 signal study. From left to right: year of observa- 
tion, effective time of good quality data used for the signal analysis, 
optical point spread function (PSF), energy threshold of the analy- 
sis, number of excess events observed and significance of the signal. 

by bad weather (including calima, i.e. Saharan sand- 
dust in the atmosphere) that limited the final data set 
and resulted in an increased energy threshold. 

All data analyzed her e were taken in th e false-source 
tracking (wobble) mode (jFomin et al.lll994l ). in which the 
telescope pointing was alternated every 20 minutes be- 
tween two sky positions at 0.4° offset from the source. 
The zenith angle of 2007 observations varied from 17° 
to 30°, in 2008 it extended up to 36°, while in 2009 it 
covered the range from 17° to 35°. 

The data w ere analyzed using the standard MA GIC 
analysis chain (jAlbert et al.ll2008at fAliu et al.ll2009f) . Se- 
vere quality cuts based on event rate after night sky 
background suppression were applied to the sample; 
28.7 hours of good quality data remained after these cuts, 
out of which 11.5 hours were taken in 2007, 8.7 hours in 
2008 and 8.5 hours in 2009. More details about the final 
data set can be found in Table [TJ For the signal study, 
a cut in the parameter size removed events with a to- 
tal charge less than 80 photo-electrons (phe) in the 2007 
data set, and 200 phe in 2008 and 2009 data sets. In the 
latter case, this cut reduces the effect of the moon light. 

Finally, for the spectrum determination an additional 
cut in PMT DC current, namely above 2.5 /LtA, was ap- 
plied to the 2009 sample , in order to reduce s ystematics 
due to the moon light (jBritzger et al.l [20 09) . resulting 
in 6.9 hours of good quality data. For the conclusive 
steps of the analysis, Monte Carlo (MC) simulations of 
7-like events were used. Hadronic background suppres- 
sion was achieved us ing the Random Forest (RF) method 
(jAlbert et al.ll200"8oT ). in which each event is assigned an 
additional parameter, the hadronness, which is related 
to the probability that the event is not 7-like. The RF 
method was also used in the energy estimation. The 
threshold of the analysis was estimated to be 80 GeV in 
2007, 150 GeV in 2008 and 160 GeV in 2009, as shown in 
Tabled 

Due to changes in the telescope performance, the sigma 
of the optical point-spread function (PSF) of 2007 and 
2008 was measured to be 13.0 mm, while in 2009 it was 
14.9 mm. To take all these differences into account, the 
data were analysed separately, using dedicated sets of 
simulated data. 

3. VHE 7 RAY RESULTS 

The 28.7 hours of good quality observations of 
PG 1553+113 carried out between 2007 and 2009 re- 
sulted i n a signal of 8. 8 a of significance according to 
eq. 17 of lLi &: Mai (j!983l ), obtained by combining the re- 
sults from each year, listed in Table [2] The signal was 
extracted by analyzing the distribution of the parameter 
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TABLE 3 

PG 1553+113 MEASURED SPECTRA 



Year 


F > 150 GeV 
[ cm~ 2 s _1 ] 


F > 150 GeV 
[Crab %] 


fo 

[ cm- 2 s- 1 TcV" 1 ] 


r 


2007 


(1.40 ± 0.38) • 10- 11 


4% 


(1.1 ± 0.3) ■ 10" 10 


4.1 ± 0.3 


2008 


(3.70 ± 0.47) • 10- 11 


11% 


(2.6 ± 0.3) ■ 10" 10 


4.3 ± 0.4 


2009 


(1.63 ± 0.45) • 10- 11 


5% 


(1.3 ± 0.2) ■ 10" 10 


3.6 ± 0.5 



Note. — Spectra of the individual years of observations of PG 1553+113. From left to right: Year of MAGIC observations; Effective 
time in hours; Integral flux above 150 GeV in units of cm -2 s _1 and Crab Nebula %, normalization factor fo in units of cm -2 s _1 TeV - 1 ; 
in the last column, the T index obtained by fitting the observed differential spectrum with a power law. The errors reported are statistical 
only. The systematic uncertainty is estimated to be 35% in the flux level and 0.2 in the power index. 



alpha, related to the incoming direction of the primary 
cosmic ray inducing the atmospheric shower. More de- 
tails on the sig nal extraction with t he alpha technique 
can be found in lAlbert et alj I2008U For the signal de- 
tection, no cut in energy was applied. 

The significance of the signal was 5.8 a in 2007, 8.1 cj 
in 2008 and 4.2 a in 2009. Due to a large difference in 
the energy thresholds and changes in the experimental 
conditions, the obtained fluxes cannot be compared di- 
rectly. A detailed spectral analysis is necessary in order 
to study the source emission. 

3.1. Integral Flux 

In order to explore the VHE 7 ray emission of 
PG 1553+113 from each year, we compared the inte- 
gral flux above 150 GeV. This value is a safe compromise 
taking into account the different energy thresholds. The 
final samples and the results of the spectral analyses are 
shown in Table |3] 

The integral fluxes measured above 150 GeV lie in the 
range of 4% to 11% of the Crab Nebula flux measured 
by MAGIC (| Albert et all l2008afl : the highest flux level 
is recorded in 2008 (0.11 Crab unitf]), a factor between 
two to three larger compared to the one measured in 2007 
(0.04 Crab units) and 2009 (0.05 Crab units). A constant 
fit to the data has a probability smaller than 0.001%. 
Such changes in the flux level observed in PG 1553+113 
are quite moderate in comparison to other monitored 
TeV blazars. For example, in Mkn 421 a flux variations 
exceeding one order of magnitude have been observed 
(je.g. Fossati et alll2008D . 

A detailed study about possible flux level variations on 
short timescale was carried out with the limiting condi- 
tion that the signal is not strong enough to allow for a 
detailed sampling on sub-day timescale. The upper panel 
of Figure [1] displays the light curve of PG 1553+113 mea- 
sured from 2007 to 2009 by MAGIC with a variable bin- 
ning. For comparison, the daily flux levels measured in 
2005 and 20 06 are shown, as extr apolated from the pub- 
lished data (j Albert et al.l l2007al) . and rescaled accord- 
ing to the power laws that interpolate the differential 
fluxes. Furthermore, the 2006 mean integral flux above 
150 GeV tak en during the multiw avelength campaign and 
reported in lAlbert eta l. (2009) is shown. The former 
data have not been used for the integral flux study, due 
to very large uncertainties related to the extrapolation 
procedure. We set 2-days, daily and monthly binning 

1 The Crab unit used in this work is an arbitrary unit obtained by 
dividing the integral energy flux measured above a certain thresh- 
old by the Crab Nebula flux m easured above the same threshold 
by MAGIC (Albert et al. 20081). 



for the 2007, 2008, and 2009 data sets respectively, ac- 
cording to the significance of the signal. The 2008 data 
are consistent with the hypothesis of constant flux with 
a probability of 50% (Figure lower panel). 

In 2007, the observed time series is consistent with con- 
stant flux (93% of probability). Nothing similar can be 
concluded for the 2009 observations since the significance 
of the signal is too low. 

In general, the high energy threshold of the analysis 
together with the weakness of the PG 1553+113 signal 
and its very steep spectrum make any variability study 
at short timescale difficult and might have hidden the 
detection of an increased activity on very short timescale. 

3.2. Differential Flux 

The differential spectra observed from PG 1553+113 
by MAGIC each year from 2007 to 2009 are shown in 
the left plot of Figure [H 

As for other blazars, each spectrum can be well fitted 
with a power law function of the form 



dF 

= f n * 

dE J0 



E 



200 GeV 



(1) 



where fo is the flux at 200 GeV and T is the power law 
index. The resulting indices are listed in the last col- 
umn of Table [3] The systematic uncertainty is estimated 
to be 35% in the flu x level and 0.2 in the power index 
(| Albert et al.ll2008al) . and is the sum of many contribu- 
tions, mainly related to the use of MC simulations in- 
stead of test beams. Thanks to the low energy thresh- 
old of the analysis of 2007 data, the corresponding spec- 
trum has a measured point below 100 GeV. In particu- 
lar, the measured differential energy flux at 98 GeV is 
(2.7 ± 0.3) TO -9 cm _2 s _1 TeV _1 , in agreement, within 
the errors, with the low energy point measured in 2005 
and 2006, (4.1 ± 1.2) • 10~ 9 cnT^s^TeV- 1 at 97 GeV. 

The 2008 differential energy spectrum measured above 
150 GeV has a slope of 4.3 ± 0.4, while the slope of 
the spectrum determined with a partially simultaneous 
sample taken during a multiwavelength campaign with 
other instruments i s 3.4 ±0.1 between 70 to 350 GeV 
(jAleksic et al.l 120101) . The different energy range char- 
acterizing the measurements fully accounts for this ap- 
parent disagreement: the spectral points measured in the 
range 150-350 GeV are, in fact, in very good agreement. 

Finally, the 2009 differential spectrum is barely deter- 
mined due to the limited signal. Except for the latter 
sample, whose significance is rather low and correspond- 
ing errors noticeably large, the power law indices describ- 
ing the spectra are compatible. This indicates that the 
shape of the emitted spectrum does not change, even if 
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Fig. 1. — Multiwavelength light curve of PG 1553+113 from 2005 to late 2009. From upper to lower panel: VHE 7 rays above 150 GeV 
measured by MAGIC (filled circles), optical flux in the R-band (triangles), X rays 0.3-10 keV flux (squares) and 7 rays above 1 GeV (open 
circles). The error bars reported have 1 a significance. Downward triangles, in lower panel, refer to 2sigma upper limits on the source flux 
above 1 GeV. 

correct for the effects due to the finite en e rgy res olution 
(procedure called unfolding, lAlbert et al.l (|2007bO 'l. The 
good agreement among these mean determinations sug- 
gests that despite the (small) variability seen on yearly 
scale, the mean flux emitted by this source is stable. A 
power law fit gives the values T — 4. 27 ± 0.14 for the in- 
dex and f = (1.61 ± 0.14) • 10~ 10 s" 1 cur 2 TcV" 1 for 
the normalization factor, with a x 2 /dof = 5.7/9 and cor- 
responding probability of 77%. The integral flux above 
150 GeV is at the level of 8% of the Crab Nebula flux. 

VHE photons from cosmological distances are ab- 
sorbed in the interaction with the EBL. Taking into 
account the EBL absorption assuming the background 
model proposed in IDom ingucz et al.l (|2011[ ) and a red- 
shift z = 0.40, the intrinsic spectrum is compatible with 
a power law of index 3.09 ± 0.20, as drawn in Figure |3l 
If we assume z = 0.45, the corresponding spectrum is 
compatible with a power law of index 2.91 ± 0.21. 
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Fig. 2. — Zoom of Figure Q] around 2008 MAGIC observations 
(lower panel). In the upper panel the corresponding optical flux 
has been superimposed for comparison. 

the total flux shows hints of (small amplitude) variability, 
as also noted for ot her BL Lacs such as 1ES 1218+304 
(jAcciari et al.ll2010D . 

The right plot of Figure [3] shows the combined dif- 
ferential spectrum of PG 1553+113 from 2007 to 2009, 
superimposed to the 200 5-2006 spectrum m easured by 
MAGIC (r = 4.21 ± 0.25. IAlbert et~aT1l2007ai) . The gray 
band represents the systematic effect on the combined 
spectrum resulting from the use of different methods to 



4. PG 1553+113 AS SEEN AT OTHER WAVELENGTHS 

Figure[T]displays the light curve of PG 1553+113 in dif- 
ferent wavelengths. The MAGIC data shown cover five 
cycles of TeV observations at energies above 150 GeV. 
The time bins used are variable, as described in the pre- 
vious section. 

The simultaneous optical R-band data are outlined in 
the second panel. These data are collected on a nightly 
basis by the Tuorla Observatory Blazar Monitoring Pro- 
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Fig. 3. — Differential energy spectra from PG 1553+113. Left Figure: comparison between 2007, 2008 and 2009 spectra. Right Figure: 
supcrimposition of 2005 -2006 spectrum, fro m Albert et al. (2007a), to 2007-2009 mean spectrum and corresponding dcabsorpti on for z = 0.4 
using the EBL model of lDo mingucz ct al. (2011). In both figures, the fit of the Crab Nebula spectrum measured by MAGIC (lAlbert et al.l 
2008a) is superimposed for comparison. 

Fermi data presented in this paper are restricted to the 
1 GeV - 100 GeV energy range and were collected from 
MJD 54682 (2008 August 4) to MJD 55200 (2010 Jan- 
uary 4) in survey mode. 

An unbinned analysis was performed to produce the 
light curve with the standard analysis tool gtlike, in- 
cluded in the Science Tools software package (version 
v09r21p00). P6_V11 DIFFUSE Instrument Response 
Functions (IRFs) were used, which are a refinement to 
previous analyses reflecting improved understanding of 
the point spread function and effective area (Abdo et al. 
2011, in preparation). For this analysis, only photons 
belonging to the Diffuse class and located in a circular 
Region Of Interest (ROI) of 10° radius, centered at the 
position of PG 1553+113, were selected. In addition, we 
excluded photons arriving from zenith angles > 105° to 
limit contamination from Earth limb 7 rays, and photons 
with rocking angle > 52° to avoid time intervals during 
which Earth entered the LAT FoV. 

A separate analysis of the high energy emission in 
each time bin was performed. All point sources in the 
1FGL within 15 degrees of PG 1553+113, including the 
source of interest itself, were considered in the analysis. 
Those within the 10 degree radius ROI were fitted with 
a power law with spectral indices frozen to the values 
obtained from the likelihood analysis of the full data set, 
while those beyond 10 degree radius ROI had their values 
frozen to those found in 1FGL. 

Upper limits at 2 a confidence level (downward trian- 
gles in Figures 1 and 4) were computed for time bins 
with Test Statistics (TS)Q < 4, and were handled as in 
the first Fermi/LAT catalog paper. The estimated sys- 
tematic uncertainty on the flux is 10% at 100 MeV, 5% 
at 500 MeV, and 20% at 10 GeV. 

As already noted, during the five years of monitoring 
the source generally showed a marginal activity in the 
VHE 7 ray band. The same behavior is followed by the 
optical flux, whose variations are limited within a fac- 
tor of four, with a maximum flux reached in 2008 and 
a minimum value in 2009. A low emission in 2009 is 
also registered at all the other wavelengths, Figure 21 



grarrfl (jTakalo et al.l 12007ft using the KVA 35 cm tele- 
scope at La Palma and the Tuorla 1 meter telescope in 
Finland. 

For the third panel in Fig 1 we used the 14 
Swift pointed observations (Gehrels et al. 2004) of 
PG 1553+113 performed from 2005-04-20 to 2010-02-05. 
We summed the data collected on 5 and 7 July 2009 in 
order to have enough statistics to obtain a good spectral 
fit. The XRT data were processed with standard proce- 
dures (xrtpipeline vO.12.6), filtering, and screening 
criteria by using the Heasoft package (v6.11). We con- 
sider the data collected in photon counting mode, and 
thus only XRT event grades 0-12 were selected. 

Source events were extracted from a circular region 
with a radius of 20 pixels (1 pixel ~ 2.36"), while back- 
ground events were extracted from a circular region with 
radius of 50 pixels away from the source region. Some ob- 
servations showed an average count rate of > 0.5 counts 
s -1 , thus pile-up correction was required. In that case we 
extracted the source events from an annular region with 
an inner radius from 2 to 7 pixels (depending on the 
source count rate and estimated by means of the PSF 
fitting technique, see Moretti et al. 2005) and an outer 
radius of 30 pixels. We extracted background events 
within an annular region centered on the source with 
radii 70 and 120 pixels. Ancillary response files were 
generated with xrtmkarf , and account for different ex- 
traction regions, vignetting and PSF corrections. We 
used the last spectral redistribution matrices in the Cali- 
bration database maintained by HE AS ARC. All spectra 
were rebinned with a minimum of 20 counts per energy 
bin to allow \ 2 fitting within XSPEC (vl2.7.0; Arnaud 
1996). We fit the spectrum with an absorbed (model 
tbabs in Xspec) log parabola law (see e.g. Tramacere et 
al. 2007), with a neutral hydrogen column fixed to its 
Galactic value (3.65xl0 20 cm- 2 ; Kalberla et al. 2005). 
The observed 0.3-10 keV fluxes obtained by Swift /XRT 
in the different observations are reported in the third 
panel of Fig. 1. 

In the lower panel, the Fermi/LAT light curve of 
PG 1553+113, computed in 10-day bins, is displayed. 

2 More information at http://uscrs.utu.fi/kani/lm/ 



3 TS is 2 times the difference of th e log(likelihood) with and 
without the source IjMattox et al.lll996T ). 
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Fig. 4. — Zoom of Figure [T] around the Cycle V observations carried out in 2009. From upper to lower panel: VHE 7 rays above 150 GeV 
measured by MAGIC (filled circles), optical flux in the R-band (triangles), X rays 0.3-10 keV flux (squares) and soft 7 rays above 1 GeV 
(open circles). Downward triangles, in lower panel, refer to 2sigma upper limits on the source flux above 1 GeV. 

suggesting that the source entered in a low activity state 
during that year, with a minimum reached few days after 
MAGIC observations. 

Figure [5] shows the result of a correlation study be- 
tween optical and TeV simultaneous observations. The 
VHE 7 ray flux above 150 GeV is plotted as a function of 
the optical flux. In order to increase statistics, we used 
for 2007/8/9 samples the daily light curve values; how- 
ever, since the optical measurements have a different time 
coverage, in some cases we derived the mean VHE flux 
from two or more consecu tive days. 2005 and 2006 data, 
from I Albert et all ()2007al) . were rejected from this study, 
due to the large uncertainty on the extrapolated flux in 
the VHE band. The mean flux yalue from 2006 mu lti- 
wavelength campaign, reported in lAlbert et al.l (|2009() . is 
included. A linear relation among the two components 
has a 74% probability, which suggests a correlation be- 
tween these two extreme energy bands. This result is 
in good agreement with the SSC model, which predicts 
a correlation between the synchrotron and the IC emis- 
sion, related to the same electron population. Due to 
the poor simultaneity of VHE data with the other wave- 
lengths, the same study has not been performed in X-rays 
and soft 7 rays. 

The X-ray light curve shows a pronounced variability, 
in contrast to optical and very high energy bands. The 
X-ray flux spans an interval of about one order of mag- 
nitude (with maximum in 2005 and minimum in 2009), 




F [mjyj 



Fig. 5. — Correlation study between PG 1553+113 optical R- 
band flux and VHE 7 ray integral flux above 150 GeV observed 
from 2006 to 2009. 

larger than that observed in the TeV, optical and GeV 
bands. The different variability displayed by the syn- 
chrotron (X-ray) and inverse Compton (GeV- TeV) com- 
ponents seems to be somewhat in contrast with the typ- 
ical behavior observed in TeV BL Lacs, showing, in ge n- 
eral, a coordinated variability (je.g. Fossati et all l^OOSfl. 
However, the sparse sampling of the observations and 

4 Rare exceptions to this rule ar e the so called "orphan" TeV 
flares, e.g. Krawczvnski et al. (2004) 
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the lack of a truly simultaneous monitoring prevents any 
strong conclusion. In particular, no optical nor gamma- 
ray data are available during the period of the maximum 
X-ray flux in 2005. Coordinated multifrequency moni- 
toring is necessary to further investigate this important 
issue. 

In a dedicated paper (|Abdo et al.l I2010blh the 
Fermi /L AT collaboration reported the analysis of the 
first year of PG 1553+113 data. They found that dur- 
ing the monitoring period the emission above 200 MeV 
is almost steady. This is in contrast with the be- 
havior of the source at higher energies (> lGeV). 
In our analysis of 2008 and 2009 LAT data (drawn 
in lower panel of Figure [T]), in fact, a steady emis- 
sion above 1 GeV has a probability smaller than 0.1% 
and is ruled out. The lowest flux, observed in April 
2009, has a value (0.5 ± 0.3) • 10~ 8 cm^s" 1 , while 
the highest flux, detected in August 2008, has a value 
(2.8 ± 0.6) • 10~ 8 cm _2 s _1 , more than 5 times higher. 
In our case, we are looking only at the upper edge of 
LAT band, probably close to the IC peak, while t he inte- 
gral flux above 200 MeV reported bv lAbdo etaLl (|2010b[ ) 
is dominated by lower energies, due to larger statis- 
tics. Therefore, we conclude that while at low energies 
(200 MeV-1 GeV) the IC continuum shows only marginal 
variability, this is not the case in the vicinity of the peak 
(some GeVs). In fact, small v ariability at GeV e nergies 
is a common feature of HBLs (jAbdo et al.H2010al) . 

5. MODELING THE SED 

In Figure we assembled the SED of PG 1553+113 
using historical data and the MAGIC spectra described 
above. Open black squares displaying radio-optical 
data are from NEE0. In the optical band, we also 
show (red diamonds) the KVA minimum and maximum 
flux measured in the period covered by MAGIC 2005- 
2009 observations together with optical -UV fluxes from 
Swift /WOT (filled black triangles, from lTavecchio et al.l 
120101) . For the X-ray data, two Swift/XRT spectra taken 
in 2005 (high flux state, red crosses, and interme diate 
state, black asterisks, from iTavecchio et al.l I2010D are 
given, and a Suza ku spectrum taken in 2006 (continu- 
ous red line, from lReimer et ai] 120081 ) . In addition, the 
average 15-150 keV flux measur ed by Swift/BAT dur- 
ing the first 54 months of survey (|Cusumano et al.ll2010f) 
is shown (black star), and the average RXTE/ ASM flux 
between March 1 and May 31, 2008 (small black square), 
from quick-look results provided by the RXTE/ ASM 
teanfl 

The green triangles correspond to the LAT spectrum 
avera ged over ~ 200 days (2008 August-2009 February) 
from lAbdo et al.l (|2010bf ) . As discussed in that paper, 
the flux above 200 MeV is rather stable, showing very 
small variability over the entire period of LAT observa- 
tions. It is likely that the variability observed at the 
highest energies is not important in determining the av- 
eraged spectrum due to the limited statistics. 

For MAGIC, we report the 2005-2006 and 2007-2009 
observed spectra (filled circles) and the same spectra cor- 
rected for the absorptio n by the EBL using the model of 
iDominguez et al.l (|2011l ) (red open circles). 

5 http:/ /nedwww. ipac.caltech.edu/ 

6 http://xte.mit.edu/asmlc/ 



We model the SED with the one-zone S SC model fully 
described inlMaraschi fc Tavecchiol ([2003). The emission 
zone is supposed to be spherical with radius R, in motion 
with bulk Lorentz factor T at an angle 9 with respect to 
the line of sight. Special relativistic effects are described 
by the relativistic Doppler factor, 5 = [r(l — /3cos#)] -1 . 
The energy distribution of the relativistic emitting elec- 
trons is described by a smoothed broken power law func- 
tion, with limits 7 m ; n and 7 ma x and break at 7b- To cal- 
culate the SSC emission, we use the full Klein-Nishina 
cross section. 

Given the large variations of the X-ray synchrotron 
flux, we decided to use the average level of the syn- 
chrotron bump as measured by XRT, including also ASM 
and BAT fluxes to constrain the model. The correspond- 
ing input parameters are listed in TablelH We also report 
the derived powers carried by the different components, 
relativistic electrons, P e , magnetic field, Pg, and pro- 
tons, P p , (assuming a composition of one cold proton 
per relativistic electron) and the total radiative luminos- 
ity L r ~ L ohs /S 2 . 

In order to investigate the role of different parameters 
in the model, we have explored their variation as a func- 
tion of the intensity of the synchrotron peak. To do so, 
we have modeled the SED considering the two extreme 
states of the synchrotron peak described above, respec- 
tively. For the SSC peak, instead, we have fixed the 
VHE data. The two curves representing the models are 
superimposed in light gray to Figure [6] The parameters 
obtained, listed in the last two columns of Table 21 are 
quite similar to the ones obtained when considering the 
average level of the synchrotron bump, except for the two 
variables B and K, the magnetic field and the electron 
density. Indeed, these two parameters regulate the rel- 
ative importance of synchrotron and SSC components. 
The state characterized by a low synchrotron emission 
has larger B and smaller K values with respect to the 
mean state modeled above. Conversely, the high syn- 
chrotron emission state has smaller values of B and a 
larger K. 

Finally, a comparison with the SED mod el obtained in 
the m ultiwavelength campaign reported in lAleksic et al.l 
(2010), reveals that the parameters used for building the 
two models are quite similar. The major differences are 
the value of the Doppler factor, which in our model is rel- 
atively higher (8 — 35) than in the previous one (5 = 23), 
and that of the magnetic field (0.5 G instead of 0.7 G). 
This difference is mainly due to the higher SSC peak fre- 
quency that we find in our data, better defined by the 
combined LAT and MAGIC spectra. 

The derived value of the total jet power, Fj Ct — 
P e + Pb + Pp = 4 x 10 44 erg/s, is consistent with 
the typical values inferre d modelling similar sources 
(|e.g. Ghisellini et al.ll20Tl . We use a relatively large 
minimum electron Lorentz factor 7 m ; n ~ 10 3 in or- 
der to reproduce the hard MeV-GeV continuum tracked 
by LAT (photon index V = 1.68 ± 0.03). The 
high value of 7 m i n implies that, as commonly derived 
in TeV BL Lacs, the relativistic electrons (and the mag- 
netic field, almost in equipartition) carry more power 
than the cold proton component. Another characteris- 
tic that PG 1553+113 shares with the other TeV BL 
Lacs is that the total luminosity L T is larger than the 
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TABLE 4 

Input model parameters for the models shown in Fig. [6] 



J. (XI ex 1 1 1 1:' L v^A 




Valup 


Value ma:E 


Value m j n 


7min 


[10 J 


2.5 


1 


5 


7b 


M r\4l 

[10 J 


3.2 


3 


1.3 


7max 


Ti n5i 

[10 ] 


2.2 


5.2 


4 . 1 


Tl \ 




2.0 


2.0 


2.0 


n2 




4.0 


3.75 


3.55 


B 


[G] 


0.5 


0.8 


0.2 


K 


[10 3 cm" 3 ] 


5.35 


3.8 


25 


R 


[10 16 cm] 


1 


1 


1 


5 




35 


35 


35 


Pe 


[10 44 erg/s] 


2.2 






P S 


[10 44 erg/s] 


1.5 






Pp 


[10 44 erg/s] 


0.34 






L T 


[10 44 erg/s] 


6.3 







a We list for the three different models plotted in Fig. [6] the min- 
imum, break and maximum Lorentz factors and the low and high 
energy slope of the electron energy distribution, the magnetic field 
intensity, the electron density, the radius of the emitting region and 
its Dopplcr factor. For the average model we also give the derived 
power carried by electrons, magnetic field, protons (assuming one 
cold proton per emitting relativistic electron) and the total radia- 
tive luminosity. 

power supplied by electrons, magnet i c field and protons. 
As discussed in lCelotti k, Ghisellini I ()2008f ). this implies 
that either only a small fraction of leptons is acceler- 
ated at relativistic energies (leaving a reservoir of cold 
pairs and/or protons) or that the jet is dissipating a 
large fraction of its power as radiation, eventually lead- 
ing to the decel eration of the flow , as in fact observed 
at VLBI scales ()Piner et al.l I2010D and envisaged in the 
models of structured jets (jGeorganopoulos fc Kazanasl 
[20031: IGhisellini et al.ll2005r i. 

6. CONCLUSIONS 

In this paper, we presented the analysis of three years 
of VHE 7 ray data of PG 1553+113 collected by MAGIC 
from 2007 to 2009. The data set was divided into indi- 
vidual years, and a significant signal was found in every 
sample, confirming PG 1553+113 as a stable presence 
in the VHE sky. The overall flux above 150 GeV from 
2007 to 2009 shows only a modest variability on yearly 
time-scale, within a factor 3, corresponding to a varia- 
tion between 4% to 11% of the Crab Nebula flux. No 
clear variability on smaller time scales is evident in the 
sample. 

For the spectral analysis, the data set was combined 
with previous observations carried out by MAGIC during 
the first two Cycles of operations, from 2005 to 2006, for 
a total of five years of monitoring. This sample was ex- 
cluded from the temporal study due to very large system- 
atics related to the flux extrapolation procedure. Despite 
the hints of variability on the flux level, the differential 
flux from each year is in very good agreement with a 
power law of constant index 4.27 ± 0.14. This behavior 
has been already obs erved in other blaz ars, such as the 
HBL 1ES 1218+304 (jAcciari et al.ll2010h . 

PG 1553+113 was also monitored in optical, X-ray and 
soft 7 ray frequencies, but only the former data could be 
used for correlation studies thanks to the large timing 
coverage. Interestingly, a hint of correlation with prob- 
ability of 74% was found between MAGIC and R-band 
optical flux levels, which in turn shows only a modest 



variability within a factor 4. A clear variability is seen 
in the X-rays and 7 rays above 1 GeV. The latter out- 
come, exploring the energies close to the IC peak, is only 
apparently in contradiction with previous results stat- 
ing a quite stable spectru m for this source i n the soft 
(> 200 MeV) 7-ray band (|Abdo et al.ll20T0bf) . The dif- 
ferent energy thresholds used in the two studies can, in 
fact, explain very well the discrepancy, as discussed in 
the paper. 

Finally, for the study of the spectral energy distri- 
bution, the mean differential spectrum measured by 
MAGIC was combined with historical data at other wave- 
lengths. Due to the large variations observed in X-rays 
and characterizing the synchrotron peak, we decided to 
use for the SED modeling the high energy bump, and 
the average level of the low energy bump. A more pre- 
cise model requires coupling the VHE 7 ray part of the 
spectrum with simultaneous coverage of the synchrotron 
peak, in particular at optical- X-ray energies. An interest- 
ing feature of PG 1553+113 is the narrowness of the SSC 
peak derived from the LAT and MAGIC spectra, imply- 
ing a relatively large value of the minimum Lorentz factor 
of the emitting electrons, 2.5 x 10 3 . This is also required 
by other HBLs with hard GeV spectra (e.g. Tavecchio 
et al. 2010). 

The MAGIC stereo system, with its increased sensitiv- 
ity and low energy threshold, is the suitable instrument 
to further investigate eventual daily scale TeV variabil- 
ity, as well as to provide a good differential spectrum 
determination below 100 GeV. 
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Fig. 6. — SED of PG 1553+113. Open black squares are radio-optical data from NED, red diamonds represent the KVA minimum 
and maximum flux measured in the period covered by MAGIC observations, together with optical-UV fluxes from Swift/WOT (filled 
black triangles). For the X-ray data, two Swift/XRT spectra taken in 2005 (high flux state, red crosses, and intermediate state, black 
asterisks) are given, and a Suzaku spectrum taken in 2006 (continuous red line). The average 15—150 keV flux measured by Swift/~BAT 
(black star) and the average RXTE/ ASM flux between March 1 and May 31, 2008 (small black square), are also represented. The LAT 
spectrum averaged over ~ 200 days (2008 August-2009 February) is given (green filled triangles). For MAGIC, we display the 2005-2006 
and 2007-2009 observed spectra (blue filled circles) and the same spectra corrected for the absorption by the EBL (red open circles). The 
average SED is modeled with a one-zone SSC model (continuous black line). Alternative SED are superimposed in light gray. Detailed 
references are addressed in the text. 
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